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Abstract 

This manuscript deals with the application of a stochastic model to face 
the risk of reactivation of latent tuberculosis infection in patients undergoing 
treatment with tumor necrosis factor inhibitors. Firstly, the paper reviews the 
approach proposed by R. S. Wallis, which consists in predicting the extent of 
side effects of a given drug through extremizing procedures on a simple set of 
parameters (such as the monthly rate of reactivation of a latent infection) . Sec- 
ondly, the paper develops an analytical analysis of this approach by stochastic 
modeling. The rate for emergence of reactivation of latent tuberculosis in- 
fection is investigated by means of Markov chains and an exact solution for 
its temporal evolution is obtained. The analytical solution is compared with 
Monte Carlo simulations and with experimental data, showing overall excel- 
lent agreement. The framework outlined in this paper is quite general and 
could be extremely promising in contributing further to detecting drug ther- 
apies side effects and could develop a (hopefully fruitful) bridge between the 
physician and the physicist viewpoints. 



1 Introduction 



Over the last decades, the improved understanding of the pathogenesis of chronic 
inflammatory diseases together with a major advance in biotechnology have accel- 
erated the development of highly selective biological therapies, designed to neutral- 
ize specific targets that mediate and sustain the clinical manifestations of diseases 
[HEIE]. These compounds, mainly monoclonal antibodies (mAb) and fusion pro- 
teins, introduced a breakthrough in the management of different conditions includ- 
ing inflammatory rheumatologic disorders [3|. In this context, the first setting of 
application of the biological agents was rheumatoid arthritis (RA), a chronic au- 
toimmune disease affecting approximately 1% of the adult population If the 
disease is not treated adequately, progressive deformity can lead to loss of quality 
of life and reduce average life expectancy by about a decade [J]. Studies on the 
pathogenic mechanisms of RA have revealed that tumor necrosis factor (TNF) is 
a cytokine playing a critical role in the inflammatory cascade that results in the 
irreversible joint damage typical of the disease [5]. Following these discoveries, a 
series of clinical trials in patients with RA showed the therapeutic benefit of TNF 
blockade [6j. As a consequence, five biological agents engineered to block TNF 
actions are currently available: infliximab, adalimumab, golimumab, certolizumab 
pegol (all of them mAb), and etanercept (a receptor fusion protein) [7]. While be- 
ing highly effective, TNF blockers have raised concerns about the potential for an 
increased susceptibility to infections, in particular the reactivation of latent tuber- 
culosis (TB) infection [8]. This is caused by Mycobacterium tuberculosis, which 
may be delighted by an extremely metastable latency lasting for decades. TNF is 
involved in both protection against mycobacterial infection and TB pathogenesis, 
since it is required for the formation of granulomas, which sequester Mycobacteria 
and prevent their dissemination [8j . The reactivation of latent TB infection has been 
associated with all TNF inhibitors, hence pre-initiation screening procedures have 
been recommended, which have successfully reduced the number of reported cases 

m 

R. S. Wallis and coworkers, in a series of papers [TU| [TT| [T2| [T3] belonging to the 
medical literature, reported the most extensive and detailed studies about the rela- 
tionship between the use of TNF antagonists and the onset of several granulomatous 
infections, with particular emphasis on TB. 

From data collected through the Adverse Event Reporting System of the USA 
Food and Drug Administration in the time- window 1998—2002, considering two test- 
case drugs with a different mechanism of action (i.e. infliximab and etanercept), one 

can check that granulomatous infections involve 54 over 10 5 patients treated with 
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infliximab and 28 over 10 5 patients receiving etanercept j 13 j . Therefore, the question 
is: As the latency in TB can last decades, are these infections (in patients under 
therapy) new ones or are they reactivation of previously encountered pathogens due 
to a suppressed immune system? The answer is by far not trivial as, for TB, there 
are no secure pathways to discriminate between a new infection or the raise of a 
previous one. Moreover, a clear methodology for finding latencies is still lacking. 

Furthermore, the rarity and different characteristic sizes of this infection in dif- 
ferent countries (ranging from e.g. 5 over 10 5 in Sweden up to e.g. 140 over 10 5 in 
Romania implies that data analysis and its subsequent interpretation must be 
carefully performed. 

In order to tackle the above-mentioned problem, one could rely on the molecular 
details of TNF processing signal, which has been (at least partially) elucidated 
(see e.g. [H]), and on a clear rationale in understanding features and limitations 
of infliximab and etanercept therapies, which could be achievable directly through 
molecular immunology approach. Beyond these "standard" strands, a completely 
different route can also be performed. 

Given the relative huge amount of collected data, the problem can also be con- 
sidered from a purely inferential viewpoint, by-passing the underlying molecular 
immunology know-how. According with this perspective, in Ref. [12], an abstract 
(logical) environment has been defined, where patients can occupy one of the five 
different states: (0) No infection, (1) New infection, (2) Latency, (3) Reactivation 
of a previous TB infection, (4) Post primary TB. Now, the patients starting the 
therapy (and hence belonging to the survey whose data have been collected) can 
correspond to either state (0) or (2), because, clearly, all the other states imply 
quantifiable sickness and the patient would then be treated for TB rather than RA. 
Then, at the end of the survey, a fraction of these patients will be in an illness state, 
i.e. either state (3) or (4). The transition rates between different states are assumed 
as free-parameters and, e.g. via Monte Carlo simulations, one can find the most 
probable rates as those able to reproduce, with the smallest error, the experimen- 
tal data. Interestingly, by comparing the parameters obtained by this extremizing 
procedure, one can conclude that the main difference between infliximab and etarne- 
cept therapies resides in the largely different rates of reactivation of a preexisting 
infection, which in case of patients treated with infliximab is one order of magnitude 
higher. 

Moreover, the approach is remarkable not only as far as RA drug usage is con- 
cerned, but also from a modeling point of view as it does not employ any biological 
or medical detail, that is, it is based solely on statistical inference and extremization 
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procedures on a (reasonable) ensemble of dataQ Consequently, the method could 
be, in principle, of widespread generality and this last point motivates further our 
investigation. 

Starting from the framework introduced by R. S. Wallis, we formalize the time 
evolution of the system as a stochastic process and we describe it in terms of a 
Markov chain. The latter is solved both analytically and numerically (through 
Runge-Kutta algorithm) and further confirmed through our Monte Carlo tests. In 
particular, the analytical solution allows to obtain a simple functional form for the 
number of patients affected by TB versus time and to evidence qualitative differences 
between different drugs; it also allows to detect a set of integrals of motion underlying 
hidden symmetries. 

2 The model 

In this section we formalize the model introduced in |12| . Such probabilistic 
model is based on purely clinical variables with the aim of reproducing data of TB 
onset in patients treated with TNF inhibitors, with particular attention on infliximab 
(an anti-TNF mAb) and etanercept (a soluble TNF receptor). 

The model, whose structure is depicted in Figure [T], consists in identifying a set 
of possible states^] for the patient subjected to drug treatments, and in fixing the 
likelihood for the patient to change his/her state. In mathematical terms, we need 
to fix the transition rates per unit time, where the timescale is set by the frequency 
of data collection, which here is one month. 

The clinical states available to a test-patient are taken as (see Figure 1) follows: 

: Absence of infection; 

1 : New infection (that after one month can give rise either to active TB or 
latent infection); 

2 : Latent infection; 

3 : Reactivated TB after latency; 

lr The last point is of fundamental importance because, in order to link the predicted probability 
obtained by the mathematical procedure with the observed frequency of the reported cases, one 
invokes the first postulate of probability theory, which states that for an infinite number of cases 
probability and frequency do coincide and the larger the sampling, the smaller the difference 
between them. 

2 Clearly, on large samples, sudden incidents may appear to some patients (e.g. death for other 
causes) or some others may assume both the drugs: the analysis has been previously purified from 
these cases [121 H5] , 
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4 : Active TB with onset one month after the infection. 

Moreover each patient is assumed to monthly change his/her state, following the 
corresponding transition probabilities, which constitute the model parameter set. 
Using t to label the time, these probabilities are: 

L : Probability of having a latent infection at the beginning of the observation 
(t — 0), while, obviously, (1 — L) is the probability of not having any infection 
at that moment; 

N : Monthly rate of TB infection during the observational time; 

P : Probability of a new TB infection to become active TB after one month; as 
a consequence, (1 — P) is the probability of this new infection to give rise to 
a latent infection; 

R : Monthly rate of reactivation of a latent TB infection; 




Figure 1: Symbolic representation of the Markov chain under investigation. 

It is worth noting that the available data consist in a collection of times correspond- 
ing to the onset of TB (in its active phase, namely a detectable scenario) after the 
beginning of the treatment with TNF blockers. As a consequence, the only states 
which are possible to observe are the states 3 and 4. 

Unfortunately, as widely discussed in the introduction, these states (that account 
for ill patients) are not distinguishable one respect to the other simply by looking 
at the data (hence motivating both Wallis study and our mathematical deepening), 
however, some progress can be made using stochastic extremization (whose most 
famous example is Monte Carlo simulation). The idea resembles the standard max- 
imum likelihood and it consists in finding the best values for free parameters such 
that the theoretical curves collapse over the experimental data [T5l IT6] . 
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As already outlined, following this procedure, the main result in [12] is that the 
principal difference between infliximab and etanercept treatments resides in different 
management of latent TB. The former seems to enhance reactivation of one order 
of magnitude more than the latter. 

2.1 Markov chains and master equations 

The model described in the previous section can be converted into differential 
equations coding for the temporal evolution of the probability of patient's states 
(which, we stress once more, can be compared to the corresponding fractions over a 
sample of patients given the huge collection of data). 

Being the states discrete, this can be accomplished in complete generality us- 
ing Markov chains (see Figure [2]), namely a (discrete-time) probabilistic framework 
where the probability of being in a given state (say B, Figure [2| at a given time 
t depends only on the probability distribution over all the states (connected to B, 
hence both A and C again in the toy-example reported in Figure [2| at the last time 
step t — 1 and on the transition rates linking these states (to B, hence Wb^-a and 

In fact, as we can see, it is instructive to start considering the toy Markov Chain 
with only three states (A, B and C), non-null transition rates wb^a and wq^b and 
time step At, shown in Figure [2j 

Note that, in the model, the probabilities of going from A to B and from C to 
B exist but not the opposite (wb^a — 0, wc->b — 0) hence, if the initial state is all 
concentrated in C, there will be no evolution, while if the starting point is spread 
among A and B, after enough time, the probability distribution will be peaked on 
C only (but in its -finite- temporal evolution resides our interest). 

It is then intuitive to understand that the probability p of being into the state B 
at time t + At is given by the probability of already being in B (hence plus 
the probability of arriving in B from A times the probability of being in A at the 
previous step (hence wb^aPa{^)) minus the probability of leaving B to C times the 
probability of being in B at the previous step (hence —w B< r -cPBif)')\ this concept 
can be written as follows: 

p B (t + At) = p B it) + Pa (t) w B ^ A At - p B (t) w c ^ B At . ( 1 ) 

Since the mathematics for continuous variable differential equations is much more 
handily and does not change significantly the perspective if the time step is small 
with respect to the global time window (in our case, the time step is one month and 

6 



the time window four years), we are allowed to consider the time as a continuous 
variable. This can be achieved straightforwardly starting by the previous equation 
using a limit procedure: 

p B (t + At) - p B (t) dp B (t) 

hm — = -— = PA{t)W B <-A - PB{t)W C ^B- 

At^o At dt 

The evolution for the probability p B (t) is then ruled by the following differential 
equation, namely a "Master Equation", which acts as a continuous counterpart of 
the Markov-chain in the discrete-time case: 

dp B (t) 

df = PA{t)W B ^A - PB{t)W C ^B ■ (2) 

In general, for a system that can be in one of N states, we need a N x N transition 
matrix w (where Wj^i is the rate for the transition from state i to state j) and the 
master equation takes the form 

dp- N N 

3=1 3=1 





Figure 2: Toy Markov chain: From the state A there is a potential flux of probability 
w b ^a toward the state B } hence we expect that, after a proper amount of 
time, some of the probability p will be drained from A to B. The same holds for 
the situation linking C to B. After an infinite time the probability of having the 
patient into the state C would be one, while being zero for the states A, B. 



3 Master equations for Wallis model 

Keeping in mind Figure [T] and the previous toy model, we can write down the 
system of differential equations describing the evolution of the five states earlier 
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introduced as follows: 



dpo 
dt 

dpi 
dt 

dp2 
dt 

dp3 
dt 

dp4 
dt 



with initial conditions 



Po(t = 0) 
p 2 (* = 0) 
Pi(t = 0) 



-Np (t), 

Np (t)-Pi(t), 

(l-P) Pl (t)-Rp 2 (t), 

Rp2(t), 

Ppl(t), 



1-L, 

Ps(t = 0) = p 4 (t = 0) 



(4) 



(5) 



0. 



The numbers indexing the probabilities mirror the enumeration of the previous sec- 
tion, hence p stands for the probability for a patient to be into the state and so 
on. 

Note that, as patients affected by active TB do not start RA therapy, we set 
Pi(t = 0) = ps(t — 0) = p±(t = 0) = 0. Furthermore, the parameter L tunes 
the initial amount of latent-TB patients with respect to free-TB patients, such that 
for L = all patients are healthy, while for L — 1 all patients display a latent 
TB-infection; 

The solution of these equations can be easily obtained using first order ordinary 
differential equations theory and the results assume the following forms: 



[1-L)e 



-Nt 



r (1-P)(1-L)N 
^ ^ (N-R)(l-R) 



Po(t) = 

Pl(t) = 
P2(t) = 

Ps(t) = 

4 (1-N)(1-R) 



R(l-P)(l-L) h 
(1-N)(N-R) \ L 

(1-P)(1-L)RN 



Rt , (1-P)(l-L)N t 
C ^ (1-N)(1-R) e 



,-Nt\ 



(i_P)(i-L)jy -m 

(1-N)(N-R) c ' 



e -) + 
1 - e" 4 ) , 
- e~ Nt + 1-N). 



j (1-P)(1-L)N 
^ ^ (1-R)(N-R) 



(i 



-Rf) 



(6) 



Of course, since the total amount of patients is conserved, C = P0+P1+P2+P3+P4 
is an integral of motion, that is 







d 
di 



C =>• C = p (t) + + p 2 (t) + p 3 (t) + p 4 (t) = const. 



(7) 
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Beyond Co, the system above admits another integral of motion Ci, namely 



dt 



P2(t) + P3(t) + ^— W) 



P - 1 

Ci = p 2 (0+P3(*) + — ^rf) = const. (8) 



This means that the rate of growth for patients in the latency branch (i.e. in 
states 2, 3) equals the rate of growth for remaining infected patients (i.e. in state 
4) weighted by a factor P~ x — 1, so that the smaller P and the larger the difference 
between the related rates. The knowledge of integrals of motion can be very useful 
as they allow to obtain information in a very simple way; for instance, if at any time 
P decreases, p±(t) must also decrease (or, analogously, P2(t) + P3{t) must increase) 
in order to maintain C\ constant. Given Cq and C±, other integrals of motion which 
are combination of C and Ci, can be trivially built. For example, C 2 = C — C\ 
fulfills 



dt 



Pl(t)+p (t) + pP4(t) 



C 2 = Pi(t) + po(t) + jp A (t) = const. (9) 



We underline that this kind of investigation can be accomplished only through an 
analytical study of the system. 

As discussed above, the fraction of active TB cases is given by the sum of the 
fraction of cases of direct TB after infection and of the fraction of cases with reac- 
tivated TB; namely, calling c(t) the total fraction of cases, we have: 

c(t)=p 3 (t)+p 4 (t). (10) 

In the above quantity, the time dependence appears only through three different 
exponential decay terms (e~*, e~ Nt , e~ Rt ), which vanish at infinite time, so that 
the solution becomes a constant term equal to 1, meaning that, if we wait for a 
sufficient long (possibly infinite) time, all patients become sick (despite obviously 
they can eventually die much earlier by all the reasons not accounted in their link 
with RA/TB). To deepen the temporal evolution of these probabilities at not so 
long times, it is useful to use a little bit of algebraic manipulation to isolate these 
constant and decay terms and rewrite the evolution in time of c(t), namely the 
fraction of activated TB cases, as follows: 

c(t) = 1 + he'* + k R e- Rt + k N e~ Nt , (11) 

where the three new constants hi, kn and kjq are related to the clinical parameters 



9 



by 





- N a-L) ( P r q-f)\ 

iV (l-N) \f -"-(l-R) J ' 




< k R 


_ j . N(l-P)(l-L) 
u ^ (1-R){R-N) ' 


(12) 


k N 

V 


- {1 - L) (P 1 R {1 - P) ) 

(l-N) \f 1 (R—N) J ■ 





It is worth noticing that in the limit of long time, the exponential terms are 
vanishing in such a way that c(t) converges to 1, as expected: 

lim c(t) = 1 , 

t— >OD 

limp 3 (t) = l-P(l-L), (13) 

t— >oo 

lim p 4 (t) = P(l - L) . 

t—>oo 

Before turning attention to the approximations that better highlight the scenario 
of smaller time windows, we stress the perfect agreement between our analytical 
solution (|6]), its numerical approximation via fourth-order Runge-Kutta algorithm, 
the Monte Carlo simulations^ and real data as shown in Figure 




months 

Figure 3: Comparison between analytical, numerical (Ruge-Kutta algorithm, ode45 
in Matlab), Monte Carlo (over a population of 10 7 virtual patients) solutions and 
real data. Note that real data are the continuous (black) curves, which are not 
represented in the legend for the sake of visibility. 



3 For a Monte Carlo simulation wc mean a simulation in which a set of virtual patients evolves 
in time following the Markov chain of Figure ^ giving a sample of the evolution of the fraction of 
cases during time. 
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Figure 4: Time evolution of the complete solution c(t) with the parameters fixed 
according to our best fits (see Table pj). Through the overlapping of the curves in 
panels (a) and (b) it is possible to note the mild role played by P. 



x 10 




•inf 
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Figure 5: Comparison among the exact solution c(t) of Eq. (11), the approximate 
solution c(°'{t) (where 0(N) terms are neglected) and the approximate solution of 
Eq. (14) where a higher order c^(t) is retained. For etanercept parameters also the 



linear approximation (15) is shown. 
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3.1 Useful approximations: What we learn by studying the 
model 



While we get the complete solution of the model, from a practical perspective, it 
would be useful to find proper approximations which would not significantly affect 
the behavior of c(t), but offer a much more handily interpretation of the results. 



In particular, Eq. (11) can be reduced to a simpler form if we note that the 



monthly rate of TB infection N is much smaller than all the other parameters, in 
agreement with studies on TB and with results found in [T2] and, a posteriori, in 
the current paper (see Table [Tf. 

Hence, as a first approximation step, we assume N ^1 and N/R <C 1 such that 
we can expand the solution at the first order in N and N/R, as follows: 



h 

k N 



[1 - L) ( P - R - — ^ ) N[l + N + 0(N 2 )] } 



1-R 



-L + (1 - P)(1 - L) 7V 



R(l - R) 

I 



1+ R + °{lp) 



-(1 - L) \ P[l + N + 0(N 2 )] + (1 - P)[l + N+ 0(N Z )} 



N _ fN 2 



Therefore, we can write the following approximation 



c(t) 



1 + (1 - L) 
(1-L) 



P-R- 



1-P 
l^R 

N 



Ne- l + 



-L + 



(1-P)(1-L)N 



1-R 



R 



-Rt 



2P(1 + N) + -(1 - P) 



-Nt 



(14) 



where we neglected terms of order 0(N 2 ),0(N 2 / R) and 0(N 2 /R 2 ). By using the 
parameter values found by Wallis to fit the data of patients treated with infliximab 
and etanercept, this approximation is extremely accurate as we can check in Figure 

El 

Let us now move further and focus on the exponential terms. First, we notice 
that e - * e~ Rt e~ Nt , due to the fact that 1 » E » Af. Moreover we can 
neglect the term e _i , as it decays very quickly and the factor k% is of order O(N), 
and we can expand e~ Nt as e~ Nt ~ 1 — Nt, given that Nt 1, namely iV 10~ 2 . 
As for e~ Rt , a similar approximation (e~ Rt ~ 1 — Rt) can be adopted as long as 
R < 1CT 2 . Thus, for times 1 ^Ct^C 1/R 1/N, we can write the following linear 
approximation 

c(t)&l + k R {l-Rt) + k N (l-Nt). (15) 



The above approximation is rather good for etanercept-treated patients, as shown 
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in Figure |5j On the other hand, if we consider infliximab-treated patients, this 
approximation does not fit data. In fact, the expansion for e~ m does not hold 
because now R m l§~ l month~ l and Rt is no longer small. 

Moreover it may be stressed that, in the regime N <C (1, P, R, L), considering R 
and P as the parameters that can be affected by the treatment (as justified above), 
the most relevant parameter in ruling the solution behavior is R. In fact by plotting 
the solution and varying the order of magnitude of P (Figure 4(b)| and 4(c) ) and 



R (Figure 4(a) ) it is possible to note the small effect of a coarse tune of P with 
respect to that of R. A qualitative argument in favor of this claim is that in the 
zero approximation of the solution (i.e. neglecting even terms O(N)) P does not 
appear at all. 





Best Fits 


L 


4.52 • 10" 4 


N 


2.88 • 10" 6 


R-INF 


2.12 ■ 1CT 1 


PiNF 


9.76 ■ 1CT 1 


R-ETA 


2.24 • 10" 2 


Peta 


8.03 ■ 10- 1 



Table 1: Best fit parameters obtained through least square method. All values are 
expressed in units months^ 1 . 

To summarize, our results confirm that the most important difference between 
drug therapies based on infliximab or etanercept is that the former enhances TB 
reactivation more than the latter: In agreement with Wallis results [12], we found 
Rinf ~ 10 Ret a- However, as N is very small (of order 10 -6 months -1 ) with respect 
to the experimental time-window (26 months), we can reasonably infer only the 
reactivation rate R, while the measure of P, is affected by a much larger uncertainty, 
so that the statement Pinf — Peta has to be claimed with more caution. 



4 Summary and outlooks 

In these notes we formalized, from a mathematical perspective, a recent series of 
papers by Wallis and coworkers |10l [TT1 [T2l H~3] . where they developed a stochastic 
approach to data analysis in drug testing. In particular, we report the application of 
a stochastic model to face the risk of reactivation of latent TB infection in patients 
undergoing treatment with TNF inhibitors. 

The approach introduced in [12] bypasses the underlying biological and medical 
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suggestions toward a completely different route, close to the general framework of 
stochastic data analysis. More precisely, the possible states of each patient are split- 
ted in five different classes (healthy with no infection, new infection, healthy with a 
latent infection, reactivation of a latent infection, active TB), which reasonably cover 
the whole plethora of events. Then, arbitrary transition rates between these states 
are introduced, implicitly defining a Master Equation and, through Monte Carlo 
simulations, the value of these rates are estimated as those that best approximate 
the experimental data. 

In the present paper we gave a mathematical backbone to this approach, building 
it on explicit Markov processes, whose continuous limit gave us Master Equations 
governing the evolution of the probability of belonging to each of the possible states, 
that we solved in all details. We developed a clear mathematical machinery able to 
recover Wallis findings, in excellent agreement with both experimental data and the 
immunological perspective. 

It is worth stressing that this methodology, being based on very standard procedures, 
has the advantage to hold beyond the test case of RA/TB. For instance, handling the 
complete (mathematical) solution allows to better account for reasonable approx- 
imations, tackling their control quantitatively (e.g. finding the proper timescales 
involved in the process or the integrals of motion constraining the evolution of the 
system). We hope that this test-case may shed light to future developments of this 
sideline approach in figuring out adverse events in drug therapies. 
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